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ABSTRACT 

The extragalactic background light (EBL) is the diffuse radiation with the second highest energy density in the Universe after the cosmic microwave 
background. The aim of this study is the measurement of the imprint of the EBL opacity to y-rays on the spectra of the brightest extragalactic 
sources detected with the High Energy Stereoscopic System (H.E.S.S.). The originality of the method lies in the joint fit of the EBL optical 
depth and of the intrinsic spectra of the sources, assuming intrinsic smoothness. Analysis of a total of ~ 10 5 y-ray events enables the detection 
of an EBL signature at the 8.8<r level and constitutes the first measurement of the EBL optical depth using very-high energy (E > 100 GeV) 
y-rays. The EBL flux density is constrained over almost two decades of wavelengths [0.30 p.m, 17 yum] and the peak value at 1.4 yum is derived as 
AF A = 15 ± 2 stat ± 3 sys nW irr 2 sr 1 . 

Key words. Gamma rays: galaxies - Cosmology: cosmic background radiation - Galaxies: BL Lacertae objects: individual: Mrk 421, PKS 2005- 
489, PKS 2155-304, 1ES 0229+200, H 2356-309, 1ES 1101-232, 1ES 0347-121 



1. Introduction 
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The extragalactic background light (EBL) is the second most in- 
tense diffuse radiation in the Universe, and its spectral energy 
distribution is composed of two bumps: the cosmic optical back- 
ground (COB) and the cosmic infra-red background (CIB). The 
former is mainly due to the radiation emitted by stellar nucle- 
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osynthesis in the optical (O) to near infrared (IR), while the lat- 
ter stems from UV-optical light abs orbed and re-radiated by dust 
in the IR domain (for a review, see lHauser & Dwe k 2001 ). 

Direct measurements of the EBL flux density prove to be 
difficult, mainly because foreground contamination, e.g. by the 
zodiacal light, can result in an overestimation. Strict lower lim- 
its have been derived from integrated galaxy counts (see, e . g . , 
iMadau & Pozzettill2000l: [Fazio et alJl2004l: iDole et alJl2006l for 
more details). The limits derived from direct measurement in the 
near IR domain typically are one order of magnitude above the 
lower limits from source counts. 

Strong constraints on the EBL density are derived using ex- 
tragalactic very high energy (VHE, E > 100 GeV)y-ray sources. 
VHE y-rays interact with O-IR photons via electron-positron 
pair prod uction, resulting in an attenuated flux that is detecte d 
on Earth riNikishovIl 963 lJelle\fll 9661: iGould & Schrederlll967l) . 
Assuming tha t there is no intrinsic break in the energy range of 
interest (as in IStecker et al.|[T992b and that the hardness of the 
spectrum is limited, stringent upper li mits on the EBL opac- 
ity to y-rays have b een derived (e.g. lAharonian et al.l l2006ct 
M azin & Raud 120071) . Studies exploiting Fermi-LAT measure- 
ments as templat es of the intrinsic spectra h a ve also recently 
been performed (iGeorganopoulos et al.l 120101: lOrr et al 1 1201 It 
iMever et al.l |2012|) . Current models of the EBL are in close 
agreement with these limits, and they converge on a peak value 
of the stellar component /IF,; ~ 12 nW irT 2 sr " 1 , yielding a con- 
sisten t value for the opacity to y-rays (see, e.g.. lDorrunguez et al.l 
120111) . 

The attenuation by the EBL is expected to leave a unique, 
redshift dependent and energy dependent imprint on the VHE 
spectra. While at energies above E > 5 to 10 TeV (depending on 
the redshift of the source) a sharp cut-off is expected resulting 
from the CIB, a weaker modulation should imprint the spectra 
in the energy range between ~ 100 GeV and ~ 5 — 10 TeV, 
resulting f rom the rise and fall of the fi rst peak of the EBL, 
the COB ( lAharonian et al.l 1 19991 120031) . A significant detec- 
tion of this modulation, localized in a relatively narrow en- 
ergy range, requires studying high-quality spectra, as, e.g., mea- 
sured during the strong flux outburst of PKS 2155-304 in 2006 
(lAharonian et al. 2007c!), under the assumption that the intrinsic 
spectra are smooth over the energy range being studied. 

This signature is searched for in the spectra of the brightest 
extragalactic blazars detected by H.E.S.S. with a maximum like- 
lihood method, leaving the parameters of the intrinsic spectra 
free. The originality and the strength of the technique lie in the 
joint fit of the EBL optical depth and of the intrinsic spectra of 
the sources, fully accounting for intrinsic curvature. This deriva- 
tion of the EBL optical depth with H.E.S.S. data does not rely 
on constraints on the intrinsic spectrum from assumptions about 
the acceleration mechanism and results in a measurement of the 
optical depth, compared to the upper limits derived in previous 
studies. 

The sample of blazars studied in this paper, the data anal- 
ysis, and the spectral fitting method are described in Sect. [2] In 
Sect. [3] the results are presented and the systematic uncertainties 
are discussed. Finally, the results of this analysis are compared 
with the current constraints in Sect. |4] 

2. Analysis of H.E.S.S. data 

2.1. Reduction of H.E.S.S. data 

The high energy stereoscopic system (H.E.S.S.) is an array of 
four imaging atmospheric Cherenkov telescopes located 1800 m 



Data set N y cr £" min - £, 

[TeV] 



Mrk421 (1) 




3381 


96.7 


0.95- 


-41 


Mrk421 (2) 




5548 


135 


0.95- 


-37 


Mrk421 (3) 




5156 


134 


0.95- 


-45 


PKS 2005-489 


U) 


1540 


25.3 


0.16- 


- 37 


PKS 2005-489 


(2) 


910 


28.9 


0.18- 


-25 


PKS 2155-304 


(2008) 


5279 


99.2 


0.13- 


- 19 


PKS 2155-304 


0) 


3499 


93.0 


0.13 - 


-5.7 


PKS 2155-304 


(2) 


3470 


116 


0.13 - 


-9.3 


PKS 2155-304 


(3) 


9555 


186 


0.13- 


- 14 


PKS 2155-304 


(4) 


4606 


132 


0.18 - 


-4.6 


PKS 2155-304 


(5) 


11901 


219 


0.13 - 


-5.7 


PKS 2155-304 


(6) 


6494 


166 


0.15- 


-5.7 


PKS 2155-304 


(7) 


8253 


191 


0.20- 


-7.6 


1ES 0229+200 




670 


12.6 


0.29- 


-25 


H 2356-309 




1642 


21.2 


0.11 - 


-34 


1ES 1101-232 




1268 


17.8 


0.12- 


-23 


1ES 0347-121 




604 


13.5 


0.13- 


- 11 



Table 1. Data sets on VHE blazars detected by H.E.S.S. that are 
used for this study of the EBL. For highly variable sources, the 
data are divided into smaller subsets that are indexed in column 
1 and correspond to restricted flux ranges. The photon excess, 
detection significance, and energy range of the spectra (in TeV) 
are given in columns 2, 3, and 4, respectively. 



above sea level, in the Khomas Highland, Namibia (23°16' 18"S, 
16°30'01"E). The Cherenkov light emitted by VHE-particle- 
induced showers in the atmosphere is focussed with 13 m di- 
ameter optical reflec tors onto ultra fast cameras (tBernlohr et al.l 
l200l lHinton1l2004l) . Each camera consists of 960 photomulti- 
pliers equipped with Winston cones to maximize the collection 
of light. The coincident detection of a shower with at least two 
telescopes improves th e y/hadron separation dFunk et al.ll2004l: 
lAharonian et al.ll2006al) . 

The data sets studied in this paper were selected with stan- 
d ard quality criteria (weat her and stability of the instruments as 
in lAharonian et al.ll2006al). and the main analysis w as performed 
with Model analysis dde Naurois & Rollandll2009l) . Based on a 
maximum likelihood method that compares the recorded images 
with simulated y-rays, this analysis improves the y/hadron sep- 
arat ion with respect to the s tandard Hillas analysis method (see 
e.g. lAharonian et al. 2006a|), especially for low energies. 

The lowest photo-electron thr eshold of 40 p.e. per camera 
image after cleaning (Loose Cuts, Ide Naurois & Rollandl l2009l) 
was adopted to cover the largest possible energy range. The on- 
events were taken from circular regions around the sources w ith 
a radius of 0.1 1°. The background was estimated with the c on- 
ventional reflected regions method (lAharonian et al 1 l2006al) . A 
minimum of three operating telescopes was required to derive 
the spectrum of a source, the redundancy allowing an improved 
reconstruction of the direction and energy of the y-rays. 

A cross-check was perfor med with the stand ard multi-variate 
analysis (MVA) described in lOhm et al 1 (120091) and an indepen- 
dent calibration, yielding consistent results. 

2.2. Sample of sources 

The detection of a subtle absorption feature, such as the effect 
of the EBL, relies on spectra measured with great accuracy, mo- 
tivating the study of the extensively observ ed, bright H.E.S .S. 
blazars. A cut on the detection significance dLi & Malll98l of 
lOcr yielded a sample of seven blazars: Mrk 421 , PKS 2005-489, 
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PKS 2155-304, 1ES 0229+200, H 2356-309, 1ES 1 101-232, and 
1ES 0347-121. 

Mrk 421 is the first extragalactic sourc e ever detected in 
the VHE energy domain (iPunch et al.lll992h . This highly vari- 
able BL Lac object is observ ed by H.E.S.S. at large zenith an- 
gles (lAharonian et al ] l2005al) . yielding a high energy threshold 
around 1 TeV but also, with a large effective area at higher en- 
ergies, photons up to ~ 40 TeV. Thu s, even considering the 
low redshift of the source, z = 0.031 dUlrich et al.lll975h . the 
EBL significantly impacts its observed spectrum, with an opti- 
cal depth t{E = 10 TeV,z = 0. 031) ~ 1. 

PKS 2005-489 (z = 0.071. | Falomo et alJfl987h and H 2356- 
309 (z = 0.165. bones etalJl2009h are two blazars at the 
~ 2% of the Crab nebula flux level, detected by H.E.S.S . 
since it went into operation (lAharonian et al ] l2006bl I2005bh . 
Whi le the latter does not show any sign of spectral variabil- 
ity dH.E.S.S. collaboration. Abramowski et al.l 1201 Oal) . an in- 
tensive observation campaign on the former revealed sig- 
nificant variations dH.E.S.S. collaboration. Acero et aH 1201 Ot 



IH.E.S.S. collaboration. Abramowski et al.l 201 ll) . 



Together with H 2356-309, 1ES 1101-232 has already been 
used for EBL studies. With a mea sured photon index sm aller 
than three for a redshift of 0.186 dRemillard et al.l Il989b . the 
spectrum of this source larg ely contributed to the st ringency 
of the upper-limit derived by lAharonian et al.l (l2006cl) . A ded- 
icated study published in 2007 did not reveal any significant flux 
variations over the obse rvation period between 2004 and 2005 
dAharonian et al.ll2007dl) . 

PKS 2155-304 (z = 0.116. iFalomo et all Il993l) is the 
brightest extragalactic source in the Southern sky, and it h as 



been w id ely studi ed with H.E.S.S. ( Aharonian et al. 2005 did . 



2007cL I2009blal : IH.E.S.S. collaboration. Abramowski et al l 
2010M 120121). It exhibited a s pectacular flux outburst in July 
2006 dAharonian et al.l l2007d) . with a flux so high that the 
number of detected y-rays exceeds by far the cumulated excess 
from all the other H.E.S.S. extragalactic sources. This study 
focusses on the high statistics data set from July 2006 and 
from a multi-wavelength campaign performed in 2008, where 
the low state of the sou rce was measured with high precision 
dAharonian et al.ll2009bl) . These detailed high quality spectra of 
PKS 2155-304 have not been used to set limits on the EBL so 
far and are responsible for the most stringent constraints derived 
in this paper. 

1ES 0229+200 and 1 ES 0347-121 are c haracterized 
by their redshift of 0.14 dSchachteret al.l Il993l) and 0.188 



( Woo et al. 20051). respect ively. The spectra of these sources 



( Aharonian et al.l 2007albh confirmed the EBL limits set by 
Aharonian et al. rd2006ch . and their light curves were compatible 
with constant flux. 

For each source, the redshift, excess, significance, 
and energy range of the detected y-rays are shown in 
the first columns of Table [TJ Blazars sometimes ex- 
hibit spectral changes correlated w i th the flux (e.g. 
IH.E.S.S. collaboration. Abramowski et al.l 1201 Obi) that could 
result in a scatter of the absorption feature estimates. Spectral 
variations can be particularly important compared to statistical 
fluctuations for highly significant (> 30c) sources. To minimize 
this effect, the data from PKS 2155-304 (high state), Mrk 421, 
and PKS 2005-489 were divided into several bins in flux with 
roughly the same logarithmic width and a similar number of 
y-rays, using data slices of 28 min duration (runs). This resulted 
in 7, 3, and 2 bins for the sources, respectively, which are 
ordered by increasing level of flux and are listed in brackets 
in Table [TJ The observational conditions for the various data 



sets on a single source vary from one set to another, implying 
different energy ranges. 

2.3. Spectral analysis 

The spectral analysis of the data sets described in Sect. 12.21 
was performed taking the EBL absorption e~ T ^ Ex '" > into account, 
where the optical depth depends on the EBL density n and on the 
energy E of the y-rays, emitted by a source located at a redshift 
Z. The E BL optical depth w as scaled with a normalization factor 
a, as in lAbdo et al] d2010l) . yielding a spectral model for each 



<f> z (E) = <Pl t (E) x exp(- ff x t(E, z, «)) 



(1) 



where 4>" M (E) is the intrinsic spectrum of the source, i.e. the de- 
absorbed spectrum assuming an EBL optical depth scaled by a. 

The template cho s en for the EBL density n is the model 
of iFranceschini et"al] d2008l) . hereafter FR08, which is repre- 
sentative of the current state of the art of EBL modelling and 
for which the optical depth is finely discretized in energy and 
redshiffl The EBL normalization factor a, defined in Eq. [TJ is 
thus an estimator of the ratio between the measured and template 
opacities T mea sured/TFR08- The particular choice of optical depth 
modelling only has a minor impact on the reconstruction of the 
EBL flux density, and the systematic uncertainty resulting from 
this choice is estimated in Sect. 13. 31 

The functional form of the intrinsic spectrum 4>" nl (E) as- 
sumed in this study is taken from very general considera- 
tions about the source physics. Blazars spectral energy distribu- 
tions are indeed commonly described with a l eptonic emission, 
e.g. w ith synchrotron self Compton models dBand & Grindlavl 
Il985h . In the VHE range, a smooth and concave spectrum is 
expected with the possible addition of a cut-off arising from 
the Klein-Nishina effect or a cut-off in the underlying elec- 
tron distribution. The conc urrent hadronic scenarios result in a 
smooth spectrum (see e.g. lMannheimlll993l; iBeall & Bednarekl 
aronianl2000l) that closely resembles the leptonic spec- 
tra in the VHE range. At the first order, the intrinsic spectra 
are described with the most natural functional form for a non- 
thermal emission: a power law (PWL), i.e. a linear function in 
log-log scale. To test for the presence of intrinsic curvature, 
the next order of complexity is readily achieved using the log- 
parabola (LP), which is the equivalent of the parabola in log-log 
scale. The exponential cut-off hypothesis is also tested (EPWL), 
since expected on theoretical grounds, and since the order of the 
equivalent log-log polynomial would be too high, unreasonably 
widening the parameter space. The next order of complexity is 
simply achieved by generalizing the last two models, adding 
a cut-off to the LP (ELP), and smoothing (y < 1) or sharp- 
ening (y > 1) the cut-off of the EPWL (SEPWL). The exact 
choice of the intrinsic models, which are detailed in Table [2] 
does not strongly affect the EBL measurement described here- 
after, as shown in AppendixlAl 

In the following, deviations from concavity are assumed to 
arise from the EBL absorption term, which is a reasonable as- 
sumption as long as the scenarios concurrent to the leptonic 
emission do not mimic the energy an d redshift depe ndence of 
the EBL optical depth (but see also iReime r! 120071 for other 



zi = 2 in steps of 6z = W' 3 and from E = 20 GeV to £1 = 170 TeV 
for 50 logarithmic steps. An interpolation in energy is perf ormed for the 
spectral analysis. 
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probes such as flat spectrum radio quasars). The energy depen- 
dence of the EBL ab sorption deviates from mere concavity (e.g. 
iRaue & Mazinl20lol) . and inflection points in the observed spec- 
tra, which depend on the redshift of the source, constitute the key 
imprint that is reconstructed in this study. 



Name 


Abbrev. 


Function 2 




Power law 


PWL 


ME/E y 


r 


Log parabola 


LP 




a-b\og(E/E ) 


Exponential cut- 


EPWL 


</>o(E/E y 


r exp(-£/£ cut ) 


off power law 








Exponential cut- 


ELP 


ME/EoT 


n -61og(£/£„) exp( _ £/ £ cut) 


off log parabola 








Super exponential 


SEPWL 




r exp(-(£/£ cul F) 


cut-off power law 









Table 2. Smooth functions describing the intrinsic spectra of the 
sources studied in this paper. 



To quantify the amplitude of the EBL signature on 
H.E.S.S. spectra, t he maximum likelihood method developed by 
lAbdo et ail d2010l) was adapted. Likelihood profiles were com- 
puted as a function of a for each data set and each smooth in- 
trinsic spectral model given in Table [2] with a ranging from to 
2.5. For each value of a, the models <p z (E) were fitted to the data 
with the intrinsic spectral parameters free in the minimization 
procedure. The best fit maximum likelihood XL was convertecQ 
into an equivalent^ 2 = -2 log X. allowing the goodness of the fit 
to assessed with the conventional^ 2 probability as a function of 
a. 

An unconventional procedure was set up to select a model. It 
ensures that the intrinsic curvature is fully taken into account and 
that the extrinsic curvature due to the EBL absorption is not over- 
estimated. Normally, the model with fewest parameters would 
be used unless a model with one extra parameter is statistically 
preferred. Here, the model with the highest x 2 probability was 
selected, regardless of the value of a for which this maximum is 
reached. Cases where two models had comparable maximum^- 2 
probabilities are discussed in the following. 

Figure Q] shows the likelihood profiles and the^f 2 probability 
profiles derived with the smooth intrinsic spectral models given 
in Table|2]for the data sets on 1ES 0229+200 and PKS 2155- 
304 (5). In the first case, the likelihood profile derived with the 
PWL (two parameters) model does not significantly differ from 
those obtained with a LP (three parameters) or an EPWL (three 
parameters), but the decrease in the number of degrees of free- 
dom with increasing complexity favours the PWL in term of x 2 
probability, as shown in the bottom panel. In this case, the con- 
ventional method and our approach select the same model. 

In the second case, corresponding to the fifth data set on 
PKS 2155-304, the LP and the EPWL significantly improve the 
fit compared to the PWL, in terms of maximum likelihood and 
of maximum^ 2 probability. The LP profile and the EPWL pro- 
file have a similar maximum likelihood and, equivalently, (since 
the two models have the same number of free parameters) a 
similar maximum x 2 probability. Instead of performing an ar- 
bitrary choice between the LP and the EPWL, the profiles of the 
SEPWL, which generalizes the EPWL, and of the ELP, which 



2 The reference energy E is set to the decorrelation energy of the 
spectrum. 

3 The log-likelihood is set to zero if the measured number of events 
matches the expected number of events in each bin. 



1 ES 0229+200 PKS 21 55-304 (5) 




Opacity normalization a Opacity normalization a 



Fig. 1. Top panels: Likelihood profiles, as a function of the opti- 
cal depth normalization for the different intrinsic models detailed 
in the legend. The examples of the data sets on 1ES 0229+200 
and PKS 2155-304 (5) are shown on the left and right, respec- 
tively. Bottom panels: Corresponding^ 2 probabilities as a func- 
tion of the optical depth normalization. The PWL and SEPWL 
models are the spectral models chosen to describe the spectra of 
1ES 0229+200 and the fifth data set on PKS 2155-304, respec- 
tively. 

generalize thes LP, were computed. According to our approach, 
the model with the highest maximum^ 2 probability, in this case 
the SEPWL, was then selected. 

As shown in Appendix lA.21 using the common criteria with 
a significance level of 2<x yields results in agreement with the 
unconventional method described here. The drawback of our ap- 
proach is a less significant measurement with regards to the con- 
ventional method, due to the widening of the studied parameter 
space and the consequently larger statistical uncertainties. The 
intrinsic spectral models that were selected for each data set are 
given in Table [3] The impact of the selection of the intrinsic 
model on the final result is investigated in Appendix IA.2I and 
is included in the systematics (Sect. [3~3l . 

3. Results 

3. 1 . Measurement of the EBL optical depth normalization 

For each data set, the likelihood -C(a) of the EBL optical depth 
normalization and of the intrinsic spectral model is compared 
to the hypothesis of a null EBL absorption £,(a = 0). The 
test statistic (TS), defined by the likelihood ratio test TS = 
2 log [X.(a = ao)/£,{a = 0)], is shown for each data set in Fig. [2] 
With a total y-ray excess of the order of 50,000 events, 
PKS 2155-304 makes a major contribution to the EBL measure- 
ment. A maximum TS superior to 16 is achieved for the data sets 
(1 ), (5), and (7), meaning that a null EBL optical depth is re- 
jected at the Act level by each of these data sets. An EBL optical 
depth scaled up by a factor two is rejected at the 3cr level by both 
the data set (6) on PKS 2155-304 and the one on 1ES 1 101-232. 
This constraint is not surprising since 1ES 1 101-232 was alrea dy 
the most constraining source used bv lAharonian et ai1(l2006cl) to 
derive an upper-limit on the EBL opacity. The bottom panel of 
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Data set 


Spectral model 


X w>) / do} 


Mrk 421 (J) 


ELP 


21.5/31 


ATI /ill / 1\ 

Mrk 421 (2) 


ELP 


46.8 / 30 


ATI /il l / ~>\ 

Mrk 421 (J) 


ELP 


34.8 / 28 


PKS 2005-489 (1) 


LP 


49.5 / 60 


PKS 2005-489 (2) 


LP 


31.8/46 


PKS 2155-304 (2008) 


ELP 


21.9 / 37 


PKS 2155-304 (i) 


PWL 


32.3 / 31 


PKS 2155-304 (2) 


SEPWL 


25.3 / 28 




SbrWL 


35.2 / 31 


PKS 2155-304 (4) 


omn tt 

SEPWL 


19.1 / 21 


FKb (j) 


hbrWL 


24.3 / 27 


I 1\ J ilJJ-JUt ^ 


T P 




PKS 2155-304 (7) 


SEPWL 


13.6/ 13 


1ES 0229+200 


PWL 


60.1/60 


H 2356-309 


LP 


70.2/61 


1ES 1101-232 


PWL 


62.6/69 


1ES 0347-121 


ELP 


31.7/35 



Table 3. Spectral modelling of the data sets used to derive the 
likelihood profiles. The spectral models (see Sect. I2.31 > are given 
in column 2, where the acronyms PWL, LP, EPWL, ELP, and 
SEPWL are explained in Table|2] The^- 2 for the best fit EBL op- 
tical depth normalization ao and the number of degrees of free- 
dom dof are given in column 3. 



Figf2]shows the less constraining contributions. Though less sig- 
nificant individually, their combination contributes to roughly a 
third of the total TS and enables a null EBL optical depth to be 
rejected at the ~ 5<x level. 

The total TS shown in Fig. [3] i.e. the sum of the individ- 
ual ones presented in Figf2] is maximum for ao = 1.27*2' j?, at 

VATS ~ 1 .80- above the unsealed FR08 template. The upper and 
lower standard deviations correspond to a variation of ATS = 1 
from the maximum test statistic TS = 77.3. The EBL optical 



depth template scaled up by a factor ao is preferred at the 8.8<x 
level to a null optical depth, where the Gaussian significance is 
approximated by the square root of the likelihood ratio test. 

No outlier is present in the set of individual profiles, with 
best fit values of 1.44+0.29 (0.6<x), 1.23±0.34(0.1cr), 1.97+0.48 
(1.5o-), 0.75 ± 0.42 (1.2o-), and 0.48 + 0.29 (1.6cr) for the five 
most constraining data sets PKS 2155-304 (1, 5, 7, 6) and 
1ES 1101-232, respectively, where the number in brackets in- 
dicate the deviation to the best fit value ao = 1.27. Similarly, 
the less constraining contributions do not differ by more than 
VTS max - TS(a-o) i5 1 -5cr from the best fit value. 

The total TS derived with the MVA analysis and an indepen- 
dent calibration is shown in Fig. [3] Though less significant, with 
a maximum TS of 33.9, the best fit value 1.24+2:^ is in close 
agreement with the optical depth normalization derived with the 
Model analysis. The larger energy range covered by the latter 
(60% wider in logarithmic scale) accounts for the difference in 
maximum TS. 



3.2. Redshift dependence 

To investigate the redshift dependence of the EBL optical depth 
normalization, the data set is divided by redshift in three groups. 
For Mrk 421 and PKS 2005-489, the TS is maximum at a(z\) = 
1.6^j j, for an average redshift of z\ = 0.051. The TS of the 
eight data sets on PKS 2155-304 (z 2 = 0.116) peaks at a(z 2 ) = 
1.36 + 0.17. With the four other data sets, a maximum TS is 
obtained for afo) = 0.71^ ^, corresponding to a mean redshift 
Z 3 = 0.170. 
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0.5 1 1.5 2 
Opacity normalization a 



2.5 



Fig. 2. Test statistic as a function of the normalized EBL optical 
depth for the intrinsic spectral models described in Table [3] The 
TS profiles are sorted by contribution to the measurement and 
the top panel shows the most constraining data sets, while the 
bottom panel shows the less constraining contributions. The ver- 
tical line indicates the best fit value derived in this study. Note 
the different scale on the vertical axis in the upper and the lower 
panel. 



h 

U 




0.5 1 1.5 

Opacity normalization < 



2.5 



Fig. 3. Combined test statistic as a function of the normalized 
EBL optical depth. The results obtained with the Model analysis 
are shown with a black line and the cross-check led with the 
MVA analysis is shown with the dashed grey line. The best fit 
value and lcr statistical uncertainties are shown with the vertical 
lines. 



Fitting the decreasing trend of the EBL normalization as a 
linear function of the redshift yields x 2 1 dof = 0.41/1, which 
does not significantly improve the fit with regards to a constant 
fit x 1 /dof = 1.83/2. A likelihood ratio test prefers the linear 
fit only at the l.lcr level. Any redshift dependence of the EBL 
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normalization in the redshift range probed is therefore neglected 
in the following sections. 

Given the limited amount of data, the deviations from the 
best fit EBL normalization ckq = 1.27*2'}? can hardly be investi- 
gated at the single data set level. For the three above-mentioned 
groups of sources, the total number of measured events in each 
energy bin (N mes ) is scaled to the expected number of events 
from the intrinsic spectra (N±, a =o)- This ratio is compared in 
Fig. |4] to the best fit model for the three average redshifts of 
0.051, 0.116, and 0.17. Abrupt changes in the amplitude of the 
statistical uncertainties (e.g. around 1 TeV for the low redshift 
group: Mrk 421 / PKS 2005-489) are inherent to the grouping of 
data sets that cover different energy ranges (e.g. the data sets on 
Mrk 421 start at ~ 1 TeV). 




1 10 

Energy [TeV] 



Fig. 4. Observed number of y-rays over number of events ex- 
pected from the intrinsic spectra vs y-ray energy. The data sets 
are grouped by similar redshift and the detected and expected 
numbers of y-rays are summed in each energy bin. The best fit 
EBL absorption is represented by the solid lines for the three 
redshifts corresponding to the groups of data sets and the shaded 
areas correspond to the ±l<x best fit EBL normalization. 



3.3. Systematic uncertainty 

An extensive investigation was undertaken of the systematic un- 
certainties arising from the method. Four sources of systematic 
uncertainties on the EBL optical depth normalization were iden- 
tified: the analysis chain (background rejection, spectral anal- 
ysis), the choice of intrinsic models and of the EBL template, 
as well as the limited knowledge of the energy scale due to the 
atmosphere. These systematic uncertainties are summarized in 
Table |4] and detailed in Appendix lAl 

The total systematic is estimated as o~ S y S (ao) = 0.25 and is 
comparable to the statistical uncertainty on the normalized EBL 
optical depth a = 1.27^-}! . 



Sources of systematics 


Estimated systematics 


Analysis chain 


0.21 


Intrinsic model 


0.10 


EBL model 


0.06 


Energy scale 


0.05 


Total 


0.25 



Table 4. Sources of systematics and estimated uncertainties 
on the normalized EBL optical depth ao = 1.27*q"J| tat . A 
full discussion of the systematic uncertainties can be found in 
Appendix lAl 



4. Discussion 

The measurement of the EBL optical depth can be converted to 
an EBL flux density, but particular attention must be paid to the 
wavelength range covered. 

A y-ray of energy E* and an EBL photon of energy e* tend 
to produce an electron-positron pa ir mostly for E*e* = (2m e c 2 ) 2 
(peak of the cross section, see, e.g. Jjauch & Rohrlichl 19761) . The 
interaction can occur anywhere along the path of the y-ray from 
the source and the relation for the EBL wavelength becomes, in 
the observer frame, 

0Wl)Um) = 1-187 x(E/l TeV) x(l+z') 2 (2) 

with z' < Z, where z is the redshift of the source and where E is 
the y-ray energy in the observer frame. To derive this relation be- 
tween the EBL wavelength and the y-ray energy, the width of the 
pair-creation cross-section as a function of energy is neglected. 
Taking it into account would result in an even wider wavelength 
coverage for a given y-ray energy range. 

The detection of an EBL flux density scaled up by a factor 
oro = 1.27+jH? stat ±0.25 sys is then valid in the overlap of the data- 
set energy ranges [(1 + z) 2 E m i n , E m!a \, where the factor (1 + z) 2 
accounts for the redshift dependency in Eq.(f2}. The measure- 
ment that is derived with all data sets is shown by the filled area 
in Fig. [5] in the wavelength range [1.2,5.5] /im, where 1.2 /mi 
(resp. 5.5 /mi) is the counterpart of the low (resp. high) energy 
bound of the Mrk 421 (resp. PKS 2155-304) data sets, as shown 
in Table 

To probe a wider wavelength range and to ensure the consis- 
tency of the modelling below and above ~ 1 /im, the TSs of data 
sets with comparable energy ranges were combined. Low EBL- 
wavelengths between 0.30 and 5.5 /mi were studied with the 
combination of the 1ES 0347-121 data set and the six PKS 2155- 
304 data sets (7, 2, 4, 5, 6, 7) while the large EBL-wavelengths 
between 1.2 and 17 /mi were probed by the 1ES 1101-232, 
1ES 0229+200, PKS 2005-489, Mrk 421, H 2356-309 data sets, 
and the two PKS 2155-304 data sets (3, 2008), all described in 
Table The normalized EBL optical depth measured in the var- 
ious wavelength ranges and the corresponding EBL flux density 
are given in Table [6] 

The lcr (statistical) contours of the EBL flux density for 
these two wavelength ranges and for the combination are com- 
pared in Fig.|5]to other measurements and limits. The first peak 
of the EBL flux density, the COB, is entirely constrained by the 
low and the high energy data sets. The systematic uncertainty is 
quadratically added to the statistical uncertainty on the measure- 
ment with the full data set in the intermediate wavelength range, 
and to uncertainties on the low and high energy measurements in 
the extended ranges. The statistical uncertainties remain domi- 
nant around 10 /im. In the UV to NIR domain, the systematic un- 
certainties, which are propagated from the optical depth normal- 
ization to the flux density as a single normalization factor, make 
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Data set 


z 


F - 


- F 


^min 


^max 






[TeV] 


Ip 


m] 


Mrk421 (1) 


0.031 


0.95 


-41 


1.2 


-49 


Mrk 421 (2) 


0.031 


0.95 


-37 


1.2 


-44 


Mrk421 (3) 


0.031 


0.95 


-45 


1.2 


-53 


PKS 2005-489 (1) 


0.071 


0.16 


-37 


0.22 


-44 


PKS 2005-489 (2) 


0.071 


0.18 


-25 


0.25 


- 30 


PKS 2155-304 {2008) 


0.116 


0.13 


- 19 


0.30 


-23 


PKS 2155-304 (i) 


0.116 


0.13 


-5.7 


0.19 


-6.8 


PKS 2155-304 (2) 


0.116 


0.13 


-9.3 


0.19 


- 11 


PKS 2155-304 (3) 


0.116 


0.13 


- 14 


0.19 


- 17 


PKS 2155-304 (4) 


0.116 


0.18 


-4.6 


0.19 


-5.5 


PKS 2155-304 (5) 


0.116 


0.13 


-5.7 


0.27 


-6.8 


PKS 2155-304 (6) 


0.116 


0.15 


-5.7 


0.19 


-6.8 


PKS 2155-304 (7) 


0.116 


0.20 


-7.6 


0.22 


-9.0 


1ES 0229+200 


0.14 


0.29 


-25 


0.45 


-30 


H 2356-309 


0.165 


0.11 


-34 


0.18 


-40 


1ES 1101-232 


0.186 


0.12 


-23 


0.20 


-27 


1ES 0347-121 


0.188 


0.13 


- 11 


0.22 


- 13 



Table 5. EBL wavelength range probed by the data sets used in 
this study. The redshifts of the sources are given in column 2. 
The energy range of the spectra (in TeV) is given in column 3, 
and the EBL wavelengths probed with the subsets are given in 
column 4, where only the peak of the pair-creation cross-section 
is taken into account. 



a non-negligible contribution to the width of the contour. The de- 
tailed study of the dependence of the systematic uncertainties on 
the wavelength, based e.g. on deviations from the EBL template 
model, is beyond the scope of this paper but the comparison of 
various modelli ngs in a complementary redshift band a nd wave- 
length range by iThe Fermi-LAT Collaboration! (1201 2b supports 
our choice of template. 



10 2 



5 



10 



I I H.E.S.S. low energy 
H.E.S.S. full dataset 
I ] H.E.S.S. high energy 

H.E.S.S. contour 
(sys + stat ) 

1 J Iniii I full I hmh ' 




10 



X [ (im ] 



Fig. 5. Flux density of the extragalactic background light ver- 
sus wavelength. The lcr (statistical) contours derived for several 
energy ranges are described in the top-right legend. The sys- 
tematic uncertainty is added quadratically to the statistical one 
to derive the H.E.S.S. contour. Lower limits based on galaxy 
counts and direct measurements are respectively shown with 
empt y upward and filled d ownward pointing triangl es (extracted 
fromlGilmore et al.ll20"l2l) . The region excluded by iMever et al.l 
(I2012T with VHE spectra is represented by the dashed area. 



r measured / T FR08 


^min 

/im 


^max 


nWnr 2 


^F, t (/l max ) 

sr 1 


1 97+O.I8 
L - Al -0.15 


1.2 - 


5.5 


14 - 8 -L7 " 


4 o+ - 6 

^•"-0.5 


1 ^4+0.19 
1 ' JH -0.17 


0.30 - 


5.5 


3.1 ±0.4 - 


4 9+0.6 


1 05 +032 


1.2 - 


17 


12.2!^ - 


q 9+1.0 
■'• -0.8 



Table 6. Measured normalization of the EBL optical depth, cor- 
responding to the lcr (statistical) contours shown in Fig. [5] The 
second column indicates the wavelength range where this mea- 
surement is valid and the third column the corresponding flux 
densities. The first line corresponds to the full data set. The sec- 
ond and third lines indicate the value derived with smaller data 
sets focussed on specific energy ranges. The systematic uncer- 
tainty on the measurements listed in the first column is 0.25. 



The contours lie in between the constraints derived with 
galaxy counts and the direct measur ements. A good agr eement 
with the VHE upper limit derived bv lMever et al.1 (1201 2l) is also 
found over the wavelength range covered, with a maximum dis- 
crepancy between 1 and 2 fim smaller than the lcr level. The 
analysis performed enables a measurement of the COB peak flux 
density of AF A = 15.0+^ ± 2.8 sys nW irT 2 sr~' at 1.4 /mi, where 
the peak value and uncertainties are derived by scaling up the 
FR08 EBL flux density by a factor ceo- This value is compatible 
with the previous co nstraints on the EBL flu x density derived 
with H.E.S.S. data bv lAharonian et aT l (l2006cT) . 



5. Summary and conclusion 

The spectra of the brightest blazars detected by H.E.S.S. were in- 
vestigated for an EBL absorption signature. Assuming intrinsic 
spectral smoothness, the intrinsic spectral curvature was care- 
fully disentangled from the EBL absorption effect. The EBL 
imprint is detected at an 8.8<x level, which constitutes the first 
measurement of the EBL optical depth using VHE y-rays. The 
EBL flux density has been evaluated over almost two decades 
of wavelengths with a peak amplitude at 1.4 fim of AFj = 
15 ± 2 sys ± 3 sys nW irT 2 sr~', in between direct measurements 
and lower limits derived with galaxy counts. 

The low energy threshold achieved with the upgrade of the 
H.E.S.S. array, H.E.S.S. II, will enable the observation of the 
unabsorbed population of y-rays and improve the constraints 
on the intrinsic spectra and thus on the absorption feature. The 
trough between the COB and the CIB will be character ized by 
the Cherenkov Telescope Array (CTA, lActis et al.ll201ll) which 
will probe energies above 50 TeV. Finally, the increasing size 
of the sample of blazars detected at very high energies will im- 
prove the constraints on the redshift dependence of the EBL and 
establish a firm observational probe of the thermal history of the 
Universe. 
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Appendix A: Study of the systematics 

The systematic uncertainties on the EBL measure ment with H.E. S.S. brightest 
blazars are investigated in this appendix. Following Sinervo 1 2003), two sources 
of systematics arising from "poorly understood features of the data or analysis 
technique" (class 2) and two sources of systematics arising "from uncertainties 
in the underlying theoretical paradigm" (class 3) are identified. The main class 
2 systematic uncertainty is evaluated with Monte Carlo simulated air showers 
imaged by the detector and passing through the whole chain analysis (see, e.g., 
Aharonian et al. 2006a, and reference therein for a description of the Monte- 
Carlo simulations). The limited knowledge of the atmospheric conditions is 
accounted for with a toy model of the detector acceptance and distribution of 
events. Class 3 systematics are characterized in this study by the choice of tem- 
plate model for the EBL absorption and the selection of the best intrinsic model 
for each data set. The impact of the latter is evaluated with the data, testing ad 
hoc intrinsic models, while the former is compared with a concurrent modelling 
established by Dominauez et al. ( 201 1). 



A. 1. Analysis chain 

Monte Carlo data (see Aharonian et al. 2006a) were used to test the analysis 
chain. Four telescopes triggered events following a PWL of photon index 2 
(hardest simulated index) were randomly removed from the simulated data set to 
create an artificial EBL attenuation. The data set studied was generated for a null 
azimuth and an off-axis angle of 0.5°. The zenith angle was fixed to 18°, close 
to the average zenith angle in the H.E.S.S. sky of PKS 2155-304, which is the 
source with the most important excess of y-rays in this study (see Sect. 12. 2) . The 
EBL optical depth normalization a was then reconstructed with these samples 
of events following a spectrum <f>(E) « E~ 2 exp(-t» X r(E, z)), where t(E, z) is 
the FR08 EBL opacity and z the redshift of the source, fixed here to z = 0.1 for 
simplicity. 

The background , particularly important for the spectral fit method described 
in Piron et "all BOOl) . was fixed to a tenth of the signal - comparable to the value 
derived for the first data set on PKS 2155-304. The reconstructed EBL normal- 
ization a is shown in the top panel of Fig. lA~T1 as a function of the simulated EBL 
normalization. The close match with the identity function strongly supports the 
reliability of the method employed. 

The parameter that seems to affect the analysis chain the most is the back- 
ground estimation, crucial for the mentioned spectral fit method. Imposing a 
background equivalent to a fiftieth of the signal, two samples of simulated events 
were studied for a null zenith and respective azimuths of 0° and 180°. The 
azimuth just indexes the data sets, since all azimuth angles are equivalent for 
a null zenith angle. The corresponding reconstructed EBL normalizations are 
represented with downward and upward triangles in the top panel of Fig. IA. II 
The associated error bars represent statistical uncertainties, related to the lim- 
ited size of the Monte Carlo samples (typically 10 4 events), that must be taken 
into account when estimating the systematic uncertainty. A first (a priori naive) 
evaluation of this systematic is the average difference a [ec0 - a simu represented 
in the bottom panel, which reads 0.17 and 0.20 for each sample. A second 
evaluation is the maximum variation in the measurement A associated with a 
Gaussian statistic, yielding one standard deviation systematics A/ Vl2 (see, e.g., 
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due to the statistics is estimated to 0.14 (variance due to uncertainties), and the 
deviation caused by the systematics is 0.10. 
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Fig. A.l. Reconstruction of the EBL normalization with Monte 
Carlo simulated air showers passing through the analysis chain. 
Three samples of Monte Carlo events are represented: the first 
one (orange squares) corresponds to the observation conditions 
of PKS 2155-304, the second and third (triangles) correspond 
to a poor background estimation. These two last sets were 
used to estimate the systematic uncertainty represented with the 
grey shaded area. Top panel: Reconstructed EBL normaliza- 
tion as a function of the simulated normalization. Bottom panel: 
Residuals, defined as the difference between the reconstructed 
and simulated optical depth normalizations. 



ISinervol[2003l) of 0.19 and 0.21, respective ly. The estimate chos en is similar to 
the excess variance estimator developed by Vaughan et al. (2003) for variability. 
Assuming that the rms difference D between the simulated and reconstructed 
values is due to both statistical and systematic uncertainties, one would write 
D 2 = 'VCci'reco _ «simu) =< °"Jtat > + °"sys> where *V is the variance estimator. We 
thus define the systematic uncertainty estimate as: 



.,)- 



(A.l) 



which reads 0.15 and 0.26 for each sample. The global systematic error using 
both samples, <x sys = 0.21, is shown in the top and bottom panels of Fig. IaTI 
This systematics estimate is similar to the two mentioned before, though a bit 
larger, which suggests a possible slight overestimation. 

To ensure that a point-to-point systematic effect does not mimic the EBL 
absorption as a function of energy, a test was performed with a bright galactic 
source, the Crab Nebula, and yielded deviations to a null EBL normalization well 
below the systematic uncertainty derived for the analysis chain. 



Fig.A.2. Likelihood profiles as a function of the normalized 
EBL opacity. The profiles were normalized to unity for clar- 
ity purposes. The dotted dashed curve is derived fitting log- 
parabolic intrinsic spectra to the data sets, while the dashed curve 
is derived by fitting exponential cut-off models. The gap between 
the two profiles due to the intrinsic spectral modelling is repre- 
sented by the grey shaded area and the double arrow. 



To ensure the reliability of the measurement, three other selection criteria of 
the intrinsic model were tested. First, the model with the best^- 2 probability was 
selected (as in the main method), but the flattest likelihood profile was used in 
case of ambiguity (e.g. between a LP and an EPWL), yielding a normalization 
of 1.18 ± 0.18, preferred at the 8.9<x level to a null opacity. A second approach 
consisted in choosing the simplest model, as long as the next order was not pre- 
ferred at the 2<x level (taking the flattest profile in case of ambiguity), yielding a 
normalization of 1 .46 ± 0. 1 1 , preferred at the 14. 3<x level to a null opacity. These 
two criteria do not change the intrinsic model for the data sets on 1ES 0229+200, 
1ES 1 101-232, Mrk 421 (2), PKS 2005-489 (1 and 2), and PKS 2155-304 (1, 6, 
and 7). A final test consisted in imposing the most complex model (an ELP) on 
the other data sets, yielding a normalization of 1.29 ± 0.18, preferred at the 7.9o" 
level to a null opacity. The above-mentioned systematic uncertainty accounts for 
the slight changes induced by the selection method and the significance of the 
result remains almost unchanged. 

It is worth noting that the particular attention paid to the intrinsic curva- 
ture of the spectra all along the analysis is not superfluous. The likelihood pro- 
file obtained assuming that the spectra are described by PWLs is maximum for 
QTowciLaw = 2.01 ± 0.07. The value derived with such a basic spectral model is 
significantly above the nominal normalized EBL opacity because intrinsic cur- 
vature of the spectra mimics the absorption effect. 



A.2. Choice of the intrinsic model 

The second systematic uncertainty arises from the choice of the model for the 
intrinsic spectra. This systematic was assessed on the data by comparing the 
total likelihood profile derived with a LP for each intrinsic spectrum, on one 
hand, and with an EPWL, on the other. The con'esponding likelihoods as a 
function of the EBL normalization are shown in Fig IA.2I where the maxi- 
mum was set to unity for clarity. The comparison of third-order models such 
as ELP and SEPWL would only drown the systematic error in the statistical one. 
The two profiles were fitted with the procedure described in Sect. [5] yielding 
«Expcut-off = l- 36 ^;?' and Q'Log Parabo la = 1.12^° }|. Using the last systematic 
estimator described in the Sect. IA.1I the difference between these two values 



A.3. Energy scale and choice of the EBL model 

The atmosphere is the least understood component of an array of Cherenkov 
telescopes such as H.E.S.S. and can affect the absoiption of the Cherenkov light 
emitted by the air showers. This absoiption leads to a decrease in the number of 
photoelectrons and thus of the reconstruct ed energy of the primary y-ray. The 
typical energy shift, of the order of 10% (Bemlohr 2000), does not affect the 
slope of a PWL spectrum, which is energy-scale invariant, but impacts its nor- 
malization. Indeed, for an initial spectrum (f>(E) = <f>o(E/Eoy r , an energy shift S 
yields ameasured spectrum (p mcs (E) = (f>o[(l +S)E/EoY r = if>' (E / Eo)~ r , where 
<f>' Q = (1 +i5)~ r </>o is the measured spectral normalization. Since the spectral anal- 
ysis developed in this study relies on the EBL absoiption feature which is not an 
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energy-scale invariant spectral model, the atmosphere absorption impact on the 
measured EBL normalization is investigated. 

A toy-model of the detector and of the atmosphere effect was devel- 
oped to account for such a systematic effect. The detector acceptance ^H(E) is 
parametrized as a function that tends to the nominal acceptance value at high 
energies, as in Eq. IA.2l 



motivates the use of the 10% energy shif0 and thus leads to a systematic error 
due to Cherenkov light absorption of 0.05. 



log 10 .?!(£) = a X [1 - bexp(-c X log 10 E)] 



(A.2) 



where ^H(E) is in m , the energy E is in TeV, and a = 5.19, b = 2.32 X 10 _ ~, c = 
3.14 are derived from the fit of the simulated acceptance. The number of events 
measured in an energy band dE is then simply AN/AE = SH(E) X <p(E) X T„i, s , 
where the observation duration T B b s was fixed to impose a total number of 
events of 10 6 . Typical event distributions for PWL and EBL absorbed PWL 
spectra are shown in the inset in Fig |A.3l A logarithmic energy binning of 
Alog| £ = 0.1 is adopted and the uncertainty on the number of events in 
each energy bin is considered to be Poissonian. To model the effect of the at- 
mosphere on the EBL normalization reconstruction, energy- shifted distributions 
AN/AE = Ji(E) X 0(£ s Mft) x T ), s were fitted with an non-shifted model, i.e. 
oc 31(E) X <p(E), with E sm = (1 + <5) X E and 0(E) cc £" r exp(-Q- X r(E,z)). 
As mentioned above the effect on the index T is null because of the energy-scale 
invariance, which is not the case for the specific energy dependence of the EBL 
opacity. A toy-model distribution that was energy shifted is shown in the top 
panel of Fig |A~3l for a redshift z = 0. 1 and an injected EBL normalization a = 1, 
corresponding to FR08 EBL modelling. The residuals Alog 10 (N evcI1 t s ) to the fit 
of a non-shifted model are shown in the bottom panel. 
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This toy model of the detector was also employed to compare independent 
EBL modellings. To probe a reasona ble range of models, the lower and upper 
bounds on the EBL opacity derived by Domineuez et al. ( 20 1 1 ) were used for the 
injected spectrum, while FR08 modelling was fitted to the event distribution. The 
variation in the reconstructed normalization is estimated to be 0.06 for a redshift 
z — 0.1, The small amplitude of the systematic effects of the atmosphere and of 
the EBL modelling choice (respectively 0.05 and 0.06) justifies a posteriori the 
use of the simple framework described in this section and does not motivate a 
deeper investigation. 
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Fig.A.3. Toy-model of the energy distribution of H.E.S.S. 
events. The inset in the top panel shows the detector accep- 
tance (black line) and the expected distributions of events for 
a PWL and an EBL absorbed PWL (green and brown lines, re- 
spectively). The injected spectra are shifted in energy to model 
the absorption of Cherenkov light by the atmosphere yielding the 
distribution of events shown in the top panel with brown filled 
circles. Fitting this distribution with a non shifted model enables 
the characterization of the atmospheric impact on the EBL nor- 
malization estimated to 0.05 for an energy shift of 10%. The 
residuals of the fit are shown in the bottom panel. 



The reconstructed and injected EBL normalizations differ by less than 0.05 
for an energy shift of 10%, while the difference can go up to 0.1 1 for an energy 
shift of 25%. The standard atmospheric conditions required by the data selection 



4 iMeve r et al. (20l3) have even shown that a precision of 5% on the 
energy scale can be achieved with atmospheric Cherenkov telescopes. 
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